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OBJECTIVE 


During  the  analysis  of  experimental  acoustic  data  from  the  non-lethal  program,  a  series  of 
technical  procedures  was  compiled  to  assist  in  a  consistent  and  standard  method  of  data  inter¬ 
pretation.  The  addressed  procedures  in  signal  analysis  represent  a  main  core  of  these  methods, 
and  are  given  with  a  representative  example  with  software  programs  (FORTRAN)  to  expedite  the 
handling  of  data.  Of  specific  concern  was  the  extraction  of  material  acoustic  impedance,  and  a 
sequence  based  on  literature  is  given  with  an  example  to  facilitate  this  extraction,  for  use  in  finite 
element  programs  for  acoustic  interaction  with  material. 


DETERMINATION  OF  COMPLEX  ACOUSTIC  IMPEDANCE 

Though  many  references  discuss  the  mathematical  and  experimental  determination  of  com¬ 
plex  acoustic  impedance,  the  approach  used  by  three  references  was  selected.  Their  method  uses 
the  common  procedure  of  determining  material  impedance  by  means  of  a  single  microphone  in  an 
impedance  tube  (refs  1  through  3). 

A  common  numerical  example  is  used  for  all  three  references,  and  English  units  are  used 
[feet,  seconds].  The  first  two  methods  are  equivalent,  though  the  mathematics  is  expressed 
differently.  The  third  accounts  for  some  attenuation  of  the  signal  down  the  tube. 

FORTRAN  programs  are  provided  in  tables  1  through  3,  with  an  example  and  input  and  output 
files  for  the  listed  program.  These  methods  use  a  sequence  of  complex  quantities  and  logarithms. 
The  third  method  can  be  solved  with  a  calculator. 

Because  none  of  the  three  references  give  a  thorough  numerical  example,  a  common  one 
was  generated  for  all  three.  Though  the  American  Society  for  Testing  Materials  (ASTM)  considers 
several  experimental  variations,  only  the  more  straightforward  method  of  two  minimum  pressures  is 
developed  in  depth.  The  ASTM  result  did  not  agree  with  the  Morse  and  Swenson  value.  In  this 
comparison,  experimental  data  is  not  being  used.  Also,  the  ASTM  approach  used  accommodates 
two  minimum  pressures,  whereas  the  others  consider  only  one. 


The  common  example  gave  the  following  values  for  the  REAL  and  IMAGINARY  parts  of  the 
impedance,  where  the  REAL  part  is  related  to  the  attenuation  coefficient.  This  coefficient  is  not  the 
surface  absorption  coefficient.  Care  must  be  given  to  note  what  specific  physical  quantity  is  being 
defined  in  the  literature. 

REAL 

IMAGINARY 

Morse  1.2630 

Swenson  1 .2630 

ASTM  1.8717 

0.7499 

0.7499 

-0.1858 

REAL  =  0  (complete  absorption  into  open  space) 
REAL  =  1  (complete  absorption  at  item) 

REAL  goes  to  infinity  (complete  reflection  at  item) 
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The  speaker  (acoustic  source)  and  sample  location,  with  general  wavelength  profile,  are 
illustrated  in  figure  1 .  During  data  acquisition,  the  microphone  is  moved  from  x=XL  down  toward  the 
speaker  at  x=0. 

SPEAKER  at  X=0  X=XC  X=XB  X=XA  X=XL(face  of  sample) 

PRESSURE:  MIN  MAX  MIN 

The  system  should  follow  these  criteria  in  setup  and  operation: 

1.  Tube  interior  must  be  straight,  smooth,  with  uniform  cross-section  and  rigid 
enough  to  keep  it  from  vibrating  during  testing. 

2.  Tube  diameter  D:  For  a  given  diameter,  acoustic  frequencies  must  be  below 
the  value  f:  f  <  (0.586*C/D)  C  =sound  speed  =  1 , 1 28  fps  at  20°C. 

In  inches:  f  <  (7930./D) 


fiHz) 

DMA 

100 

79 

500 

16 

For  a  4-in.  diameter  pipe,  the  highest  frequency  is  2,000  Hz,  and  12,000  Hz  for  a 
1-in.  diameter.  However,  in  practical  use: 

Inside  diameter  (in.)  Use  frequencies  (Hz)  below 

4  1500 

1  6,000 

3.  Tube  length  L:  For  measuring  at  least  two  minimum  pressures  (if  the  ASTM 
method  is  used),  and  to  avoid  any  measurement  within  one  tube  diameter  from  the 
speaker,  the  length  should  be:  L  >  [(0.75*C/f)  +  D] 

For  a  frequency  f=25  Hz,  and  D=4  in.,  L  >  34  ft 
For  a  frequency  f=50  Hz,  and  D=4  in.,  L  >  17  ft 
For  the  other  methods,  a  shorter  tube  would  satisfy 

4.  Specimen:  The  specimen  with  its  detachable  holder  must  make  an  airtight  fit  at 
the  tube  end.  The  specimen  should  be  flat  and  of  the  same  cross-section  of  the 
tube.  The  holder  should  be  a  massive  sound-reflective  metal.  Two  or  more 
specimens  should  be  measured  and  the  results  averaged. 

5.  Speaker:  The  speaker  may  face  directly  into  the  tube  or  be  mounted  in  an 
enclosure  to  the  side  if  it  is  too  large.  Some  absorbent  material  should  be  placed 
at  the  speaker  end  to  avoid  standing  wave  resonance.  The  speaker  frequency 
should  be  determined  to  within  1%,  or  within  1  Hz. 
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6.  Microphone:  If  the  microphone  is  small  enough,  it  can  be  placed  inside  the  tube, 
and  moved  by  a  rod  connection  along  the  length  of  the  tube.  Or  the  microphone 
can  be  attached  to  a  probe  extended  radially  to  the  tube  center  from  a  longitud¬ 
inally  slotted  tube.  An  axial  position  scale  measurable  to  the  nearest  millimeter  is 
needed. 

A  common  example  to  determine  material  impedance  is  used  for  all  three  approaches. 
METHOD  A.  Vibration  and  Sound  (ref  1) 

Experimental  data  needed:  XL,XA,XB,PMIN,PMAX 

This  approach  does  not  account  for  energy  loss  by  sound  waves  through  absorption  in  the 

tube. 

X=0  =  location  of  speaker 

X=XL  =  distance  from  speaker  to  specimen  face  =  30  ft 
X=XA  =location  of  pressure  MIN  from  specimen  =  26  ft 
X=  XB  =  location  of  pressure  MAX  =  20.36  feet. 

XLAM  =  wavelength  =  22.56  ft  (XLAM=  C/f  =  1,128  fps  at  200/  50  Hz  signal) 

PMIN  =  pressure  at  XA  =  0.1  Pascal  or  volt 
PMAX  =  pressure  at  XB  =  0.2  Pascal  or  volt 
PI  =  3.141592654 

Find  ALPHA  and  BETA  from 

BETA  =  1.-2.*(XL-XA)/XLAM  =  1  .-2.*(30-26)/22.56  =  0.64539  (1) 

tanh(PI*ALPHA)  =  PMIN/PMAX  =  0.170.2  =  0.5  (2) 

This  method  is  accurate  when  ALPHA  falls  between  0.02  and  0.5. 

Before  equation  2  is  solved,  the  following  definitions  illustrate  the  relevance  to  the  real 
impedance  THETA  that  is  used  in  ANSYS  finite  element  attenuation: 

ZETA  =  Average  specific  acoustic  impedance  of  item  at  XL 
ZETA[pressure/velocity]  =  ZM[force/velocity]/S  (area  of  tube) 

ZM  =  mechanical  impedance  of  item  at  XL 
XSI  =  Dimensionless  acoustic  impedance  ratio  at  XL 

XSI  =  Average  specific  acoustic  impedance  at  XL  /fluid  characteristic  impedance 

XSI  =  ZETA/(R*C)  =  THETA  -  i*CHI  =  TANH(PrALPHA  -  i*PI*BETA)  (3) 

XSI  =  P/(U*R*C) 

THETA  =  ACOUSTIC  RESISTANCE  (dimensionless)  =  real  part  of  XSI. 

CHI  =  ACOUSTIC  REACTANCE  (dimensionless) 

R*C  =fluid  characteristic  impedance  (density  *  sound  speed) 

U  =  particle  velocity 
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ANSYS  employs  the  parameter  MU: 

MU  =  1/THETA  =R*C  (of  fluid)/Real  part  of  material  impedance  =  R*C/(R*C*THETA) 

For  MU  =  0  =  high  impedance  at  item.  Item  boundary  reflects 

MU  =  1  =  all  absorption  at  the  item  boundary 
MU  around  E8  for  open  space 

Solve  equation  2  by  a  calculator.  This  gives  PrALPHA  =  0.549306 

Solve  equation  3  by  complex  exponentials  in  Morse  FORTRAN  program  (table  1). 

XSI  =  tanh  S  =  sinh  S/cosh  S  =  (e**S-e**(-S))/(e**S  +  e**(-S)) 

where  S  =  PI*ALPHA  -  i*PI*BETA 

XSI  is  expressed  in  FORTRAN  as  complex  natural  exponents: 

XSI  =  tanh  S  =  [cexp(S)  -  cexp(-S)]/[(cexp(S)  +  cexp(-S)] 
where  XSI  and  S  are  declared  complex. 

For  our  case:  XSI  =  tanh  (0.5493  -  i*2.0275),  S  =  S(0.5493, -2.0275) 

With  this  complex  value  for  S  in  Morse  FORTRAN  program,  the  real  part  of  XSI  is  THETA  (the 
acoustic  resistance)  and  the  imaginary  part  of  XSI  is  the  negative  of  CHI  (the  acoustic  reactance). 

For  this  example:  THETA  =  1 .2630,  CHI=  0.7499 

The  absorption  coefficient  is  easy  to  determine. 

Find  standing  wave  ratio:  SWR  =  PMAX/PMIN 

SWR  =  1  =  ALL  ABSORBED  (PLANE  WAVE). 

SWR  =  INFINITY  =  NOTHING  ABSORBED. 

R  =  REFLECTION  COEFFICIENT  =  ((SWR-1)/(SWR+1))**2  =  (1.3)**2=  .111 
IF  R  =  1,  ALL  REFLECTED. 

AB  =  ABSORPTION  COEFFICIENT  =  1.-  R  =  0.888 

METHOD  B.  Impedance  from  Principles  of  Modern  Acoustics  (ref  2) 

Experimental  data  needed:  XL,XA,PMIN,PMAX 

1 .  Same  physical  setup  as  for  METHOD  A. 

2.  Find  A  from:  A  =  (XA-XL)/XLAM  where  wavelength  =  XLAM  =  C/f,  and 
C=1 128  ft/s  at  20°C,  and  atmospheric  pressure. 

3.  Usef=50Hz.  XLAM=  22.56  ft,  XL=30  ft.  XA=26  ft. 
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4.  A=  (26-30)/22.56  =  -0.1773  Minimum  amplitudes  relate  to  phase  by:  4*PI*f*(XA- 
XL)/c  =  (2*n+1)*PI-PHA.  (n=0,1,2,3...) 

PHA  =  phase  angle  of  P(reflected)/P(incident). 

First  min  is  related  to  phase  PHA  of  wave: 

PHA  =  PI  -  4*PI*F*(XA-XL)/XLAM  =  Pr(1-4*A) 

PHA  =  Pl*(1  -  4(-0.1773))  =  5.3697  radians  =  307.7  degrees 

5.  Find:  (SWR-1)/(SWR+1)  =  1/3  =  Q  (where  SWR  =  PMAX/PMIN=.2/.1=2) 

6.  Complex  impedance  =  Zo  *[(  1  +Q*exp*(i*PHA)]/[1  -Q*exp(i*PHA)] 

Complex  impedance  =  Zo*[(1+Q*exp(i*5.3697)]/[1-Q*exp(i*5.3697)] 

Zo=  R*C  =  413  [kg/(m**2-sec)]  (not  needed) 


Use  complex  exponential: 

THETA  -  i*CHI  =  Complex  impedance/Zo  =  [1+Q*expS]/[1=Q*expS],  where  complex  S 
only  has  an  imaginary  part  =(0,5.3697) 

Then  XSI  =  THETA  -  i*CHI  =  [1+Q*CEXP(S)]/[1-Q*CEXP(S)]. 

Use  the  Swenson  FORTRAN  program  to  determine  THETA. 

Result:  THETA  =  1.2630  CHI  =0.7499 

METHOD  C.  ASTM  Method  C384-90A  (one  maximum  and  two  minima) 

Data  needed:  XLAM,XA,XC,PMIN1,PMAX,PMIN2,  with  optional  calibration.  Same  conven¬ 
tion  of  x=0  at  speaker  location  is  used  as  in  METHODS  A  and  B. 

This  method  requires  at  least  two  minimum  pressures  be  detected. 

Here  the  tube  LENGTH  should  be:  L  >  0.75*XLAM  +  D.  (D=tube  diameter) 

Procedure:  (Use  XLAM(wavelength)  =  22.56  ft  for  a  50  Hz  signal). 

1.  Place  METAL  END  REFLECTOR  at  x=XL,  the  position  where  the  specimen  face 
will  be. 

2.  Measure  XAC,  the  X  coordinate  for  first  minimum. 

3.  Measure  XCC,  the  x  coordinate  for  second  minimum. 

4.  XL  calculated  as  =  (3*XAC-  XCC)/2. 

5.  Place  front  of  SPECIMEN  where  front  of  the  reflector  had  been. 

6.  Determine  coordinates,  XA  and  XC,  of  the  first  two  minimum  pressures  from 
specimen. 

7.  Calculate:  D1=  XL-XA  =  distance  of  first  minimum  to  face  of  specimen. 

8.  Calculate:  D2  =  XA-XC  (which  is  about  half  wavelength). 

9.  Adjust  signal  so  minimum  pressures  are  more  than  10  dB  above  the  noise. 

1 0.  Don’t  measure  MAX  at  face  of  specimen. 

1 1 .  PMIN1  =  first  min  pressure.  PMAX  =  first  maximum  pressure.  PMIN2  =  second 
minimum  pressure.  These  can  be  in  Pascals  or  volts  (one  or  the  other). 
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12.  Transform  these  pressures  to  dB: 

13.  PMIN1  =20.*log(PMIN1/29E-6).  PMAX=20*LOG(PMAX/29E-6) 

14.  PMIN2  =20.*log(PMIN2/29E-6). 

15.  Calculate:  L1=  PMAX-PMIN1 

16.  L2=  PMAX-PMIN2. 

17.  L0=  LI  +  (L1-L2)/2 

1 8.  Calculate  Ko=  1 0**(L0/20) 

19.  Calculate:  M={Ko  +  1/Ko)/2 

20.  Calculate:  N  =  (Ko  -  1/Ko)/2 

21 .  Reflection  phase  angle:  PHI  (radians)  =  (D1/D2  -  0.5) 

Definitions  for  r  (real  part  of  impedance),  and  x  (imaginary  part  of  impedance): 

1.  The  impedance  ratio:  z/(rho*c)  =  r/(rho*c)  +  i*x/(rho*c) 

2.  Where:  THETA  =  1/(M-N*cos  PHI)  =  r/(rho*c)  =  resistance  ratio 

3.  Where:  CHI  =  (N*  sin  PHI)/(M-N*cosPHI)  =  x/(rho*c)  =  reactance  ratio 

Notes: 

1 .  If  D1/D2  <  0.5,  phi  is  negative,  and  x/(rho*c)  is  negative. 

2.  Specific  admittance  =  1/z  =  y 

3.  Specific  normal  acoustic  admittance  ratio  =  y*rho*c 

4.  Admittance  ratio  (y*rho*c)  is  larger  than  specific  admittance  (y)  by  about  400  in  SI  units. 

Example:  No  specimen  is  present,  only  the  metal  reflector. 

XA=30-XLAM/4=  24.36 
XC=XA-XLAM/2  =  13.08 
XL=(3*XA=XC)/2  =  30  ft 

The  specimen  is  placed  in  position: 

XA=26 

XC=14.72 

D1=XL-XA=30-26=4ft 
D2=XA-XC  =11.28  ft 

Suppose  pressures  are:  PMAX,  PMIN1,  PMIN2  =  0.2,  0.1,  0.09.  Convert  to  dB  gives:  76.77, 
70.75,  69.83  =  (20*log(P(peak)/29E-6) 

LI  =PMAX  -  PMIN1  =  76.77  -  70.75  =  6.02  dB 

L2  =PMAX  -  PMIN2  =  76.77  -  69.83  =  6.94  dB 

LO  =  L1+(L1-L2)/2  =  5.56  dB 

K0  =  10**(L0/20)  =1.897 

M  =  (K0+1/K0)/2  =1.2122 

N  =  (K0-1/K0)/2  =  0.6852 

PHI  =  (D1/D2  -0.5)  =  -0.1454  RADIANS 

THETA  =  1/(M-N*cos  PHI)  =  1.8717 

CHI  =  (N*sin  PHI)/(M-N*cos  PHI)  =  -0.1858 
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DETERMINATION  OF  SINGLE  PULSE  SPECTRAL  AND  ENERGY  CONTENT 

The  following  procedure  covers  the  essential  steps  in  pulse  analysis,  with  illustrations  made 
with  Vu-Point  3  software  for  an  experimental  impulse  noise. 

1 .  Scale  the  signal  to  actual  Pascal  pressure,  and  plot  it  lineariy  in  pressure  and 
time.  Determine  the  peak  pressure. 

2.  Determine  the  RMS  value  of  the  signal.  This  is  generally  not  well  defined  for 
pulse  signals.  An  arbitrary  interval  is  from  signal  start  to  the  time  where  the 
signal  finally  stays  within  the  boundary  of  plus  and  minus  ten  percent  of  the  peak 
pressure  (20  dB  down  from  the  peak). 

3.  In  preparation  for  the  FFT  (fast  Fourier  transform),  use  a  Hanning  (A=0.50)  or  a 
Hamming  (A=0.54)  filter  on  the  N  data  points:  [A+(1-A)*cos(2*pi*n/N)].  Also 
remove  linear  trends  in  the  signal.  For  a  pulse  with  ample  time  before  and  after 
the  pulse,  use  or  non-use  of  the  filter  is  indifferent.  This  ample  time  should  be 
present  to  avoid  signal  distortion  if  the  filter  is  used.  Obtain  the  FFT  of  the  signal. 

4.  Obtain  the  real  spectral  ENERGY  DENSITY,  FF*,  which  is  the  FFT  of  the  auto¬ 
correlation  function. 

5.  Integrate  the  FF*  to  obtain  a  step  function  showing  relative  energy  content  as  a 
function  of  frequency. 

6.  If  there  are  dominant  frequencies,  find  the  percent  contribution  to  intensity 
(energy)  from  each  dominant  frequency  by  dividing  each  step  increase  by  the 
overall  maximum.  First  subtract,  if  present,  any  energy  level  around  zero 
frequency. 

This  procedure  is  illustrated  by  these  figures: 

Figure  2  is  a  time  signal  from  440  ms  to  820  ms.  The  sampling  rate  (DT)  is  20.833  \is 
(48,000  data/s).  The  recording  time  (T)  is  0.38  sec.  Number  of  data  points  (N)  is  =  T/DT  =  0.38/ 
20.833  jis  =  18,240. 

If  all  the  data  points  are  used  in  the  FFT,  then  in  the  FFT: 

Frequency  resolution  (DF)  =  1/1=  1/(N*DT)  =  2.63  Hz 
Frequency  range  (F)  =  1/(2*DT)  =  24,000  Hz 

If  the  analysis  program  uses  a  number  of  data  points  N  =2**n  (n=..12,13,14..), 
then  the  program  chooses  n  to  have  2**n  close  to  N.  As  examples: 

2**14  =  16,384 
2**15  =  32,768 
2**16  =  65,536 
2**17  =  131,072 
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The  program  may  choose  2**14  points  and  eliminate  data  from  the  18,240,  or  use 
32,768  and  pad  the  program  with  data.  The  frequency  resolution  (DF)  of  the  resulting  FFT  can  be 
checked. 


Figure  3  is  the  FFT  of  figure  2.  The  frequency  resolution  is  1 .46  Hz.  Here,  the  program 
internally  chose  2**15  data  by  adding  14,528  points.  A  short  recording  time  T  for  a  pulse  should  be 
avoided.  Dimensions  are  [Pascal/Hz]. 

Figure  4  is  the  real  spectral  energy  density  of  the  signal  as  a  function  of  frequency. 
Dimensions  are  [Pascal/Hz]**2. 

Figure  5  is  the  integral  of  the  spectral  density,  showing  accumulative  energy,  as  a 
function  of  frequency.  Dimensions  are  [Pascal]**2. 

Figure  6  is  the  INVERSE  of  the  FFT  in  figure  2.  This  fairly  well  reproduced  the  time 
signal.  Dimensions  are  [Pascals]. 

Figure  7  is  an  FFT  with  the  Hanning  filter  of  the  immediate  pulse  itself.  Though  it 
resembles  the  FFT  in  figure  3,  its  frequency  resolution  is  coarse  at  5.06  Hz.  The  recording  time  T  of 
only  the  pulse  length  is  insufficient. 

Figure  8  is  the  INVERSE  of  the  FFT  in  figure  7.  Use  of  the  filter  has  appreciably  dis¬ 
torted  the  pulse  representation. 

Figure  9  is  the  original  time  pulse  of  figure  2,  but  with  a  broader  time  base  of  2  sec.  The 
number  of  data  N  here  =  2/20.833  ps  =  96,000. 

Figure  10  is  the  FFT  with  a  Hanning  filter  of  figure  9.  The  frequency  resolution  is  0.366 
Hz.  The  program  chose  2**17  and  added  35,072  points  to  the  available  96,000  for  a  total  of 
131,072. 


Figure  11  is  the  FFT  of  figure  9,  but  no  Hanning  filter  is  used.  This  result  is  essentially 
the  same  as  figure  10,  since  the  signal  in  figure  9  is  isolated  in  time  in  both  directions. 


ACOUSTIC  SIGNAL  CONVERSIONS 

1 .  Determination  of  peak  pressure  of  discrete  signals  in  a  composite  signal: 

If  the  fast  Fourier  transform  is  not  cluttered  and  shows  a  signal  composed  of  discrete 
frequencies,  then  the  peak  pressure  in  Pascals  of  each  discrete  frequency  is  obtained 
by  multiplying  by  two  (waves  are  positive  and  negative  by  FFT  algorithms)  the 
integrated  area  over  each  frequency  occurrence  in  the  FFT.  This  procedure  is  NOT 
used  if  a  cluster  of  frequencies  occurs  around  the  various  dominant  frequencies. 
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2.  Determination  of  pressure  at  distances  from  the  source: 

Estimating  acoustic  signal  dB  at  various  distances  from  the  source  must  be  done 
carefully.  The  effect  of  humidity,  fog,  ground  hardness  or  porosity  (snow),  ground 
reflection  and  temperature  gradient  in  the  region  above  the  ground,  and  size  of  the 
emitting  source,  have  been  reported,  but  with  limited  usefulness  for  some  engineering 
applications.  The  final  recourse  is  to  do  calibrated  acoustic  testing  with  sources  and 
along  distances  where  subsequent  investigation  is  relevant. 

The  calculation  of  energy  decreasing  as  inverse  distance  squared  (pressure  decreasing 
as  inverse  distance)  from  a  point  source,  into  a  non-attenuating  sphere,  is  simple,  but 
not  realistic: 

This  assumes  intensity  varies  as  inverse  distance  squared  along  a  radial  line  between 
source  and  microphone: 

12  =  11  *  (R1/R2)**2 

11  =  intensity  at  radial  distance  R1  (microphone  location)  from  source. 

12  =  estimated  intensity  at  radial  distance  R2  from  sound  source. 

If  II  =  1.79  W/m**2  or  122.54  dB  at  R1  =  25  meters,  then  at  R2=100  meters: 

12  =  1.79  (25/100)**2 
l2  =  0.11188W/m**2. 

dB(at  distance  R2)  =  10  log  (.1 1 188/E-12)  =  1 10.49  dB 

However,  if  the  pressure  at  the  microphone  is  calibrated  at  two  or  more  places,  the 
distance-pressure  data  can  be  fit  to  a  geometrical  regression  to  obtain  coefficients  A 
and  B  to  predict  pressure  at  a  location  X  that  is  in  a  reasonable  vicinity  where  calibration 
was  done:  P=A*(X**B)) 

If  two  calibration  points  are  made:  P1=  2  Pascals  at  x1=  3ft.  and  P2=  0.3  Pascal  at 
x2=10  ft,  then  the  coefficients  are: 

B=log(P2/P1)/log(x2/x1)  =  log(.3/2)/log(10/3)  =-1.576 
A=P1/(x1**B)  =2/(3**-1.576)  =11.297 

At  a  distance  X=15  ft,  the  estimated  pressure  is  P  =  11.297*(15**-1.576)  =  0.158  Pa. 

For  three  or  more  points,  a  standard  curve  fit  is  used  to  generate  A  and  B  for  the  pres¬ 
sure  fit. 

3.  Adding  dB  from  composite  sine-type  signals: 

Four  frequencies,  each  with  a  peak  Pascal  pressure  amplitude,  are  added  to  form  a 
composite  signal: 
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Input  amplitude 


Freauencv 

Pascal 

1000 

10 

110.75 

1500 

5 

104.73 

2000 

1 

90.75 

4000 

5 

104.73 

dB  is  from:  dB=  20  log  (P/29E-6)  [P  in  peak  Pascals] 

This  composite  signal  has  a  Pascal  RMS  of  8.69.  The  signal  dB=  20  log(P[RMS]/20E- 
6)  =  20  log(8.69/20E-6)  =  112.76  dB 

If  the  peak  amplitudes  (or  dBs  of  the  COMPOSITE  signal)  are  known,  take  the  highest 
value,  and  find  the  difference  of  the  others: 

Highest  dB  =  110.75 

2nd-highest:  1 10.75-104.73  =  6.02  =DIF  DB 
3th-highest:  1 10.75  -  104.73  =  6.02  =DIF  DB 
4rd-highest:  110.75-90.75  =  20.0  =DIF  DB 

Use  the  quadratic  expression,  ADD  DB,  to  find  what  additional  dB  is  added  to  the  dB  of 
the  dominant  frequency,  by  each  frequency  other  than  the  dominant  one: 

ADD  DB  =  3.0  -  0.427*(DIF  DB)  +  0.0153*  (DIF  DB)**2 

In  this  case: 

2nd:  ADD  DB  =  3.0  -  0.427(6.02)  +  0.0153*(6.02)**2  =  1.0 
3rd:  ADD  DB  =  3.0  -  0.427(6.02)  +  0.0153*(6.02)**2  =  1.0 
4th:  ADD  DB  =  0.  If  DIF  DB  is  14  or  more,  there  is  no  contribution. 

Total  dB  of  signal:  110.75+1+1  =  112.75.  This  is  close  to  the  known  dB  of  the  entire 
signal. 


Pressure-Power-dB  Conversions: 

Pressure 

(volts) 


Power 

(watts) 


dB 

fadditive  increment) 


1.0 

1.0 

0.0 

1.122 

1.26 

1.0 

1.26 

1.59 

2.0 

1.414 

2.0 

3.0 

2.0 

4.0 

6.0 
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Pressure 

Power 

(volts) 

(watts) 

dB 

(multiolvina  increment) 

(multiolvina  increment) 

(additive  increment) 

2.828 

8.0 

9.0 

3.162 

10.0 

10.0 

4.0 

16.0 

12.0 

5.0 

25.0 

14.0 

10.0 

100.0 

20.0 

100.0 

10,000.0 

40.0 

1,000.0 

1,000,000.0 

60.0 

A  12  dB  increase  of  an  80  W  signal  is  16*80  =  1280  W. 

A 12  dB  increase  and  3  dB  decrease  is  8*80  =  640  W. 

A  twofold  increase  in  voltage  results  in  a  6  dB  signal  increase. 
Conversely,  the  increments  can  be  dividing  and  subtractive. 

5.  Intensity,  Pressure  (Peak,  RMS)  relations: 


Effect 

Intensity 

(dB) 

Intensity 

(W/m**2) 

Peak  pressure 
(Pascal) 

Peak  pressure/SQRT2 
(Pascal  RMS) 

— 

194 

25.1E6 

143,340.0 

101,360.  (atm,14.7psi) 

191 

12.6E6 

101,360.0 

71,700.0 

— 

180 

1E6 

28,900.0 

20,400.0 

~ 

160 

10,000.0 

2890.0 

2040.0 

~ 

140 

100.0 

289.0 

204.0 

Pain 

120 

1.0 

28.9 

20.4 

Riveter 

95 

0.0032 

1.631 

1.15 

Train 

90 

0.001 

0.91 

0.644 

Busy  Street 

70 

E-5 

0.091 

0.0644 

Quiet  radio 

40 

E-8 

0.0029 

0.002 

Leaves 

10 

E-11 

9.1  E-6 

6.4  E-6 

Threshold 

0 

E-12 

28.9  E-6 

20.4  E-6 

6.  Acoustic  wavelengths  for  20C  and  1  atmosphere: 

Wavelength  for  frequencies  based  on  speed=344  m/s  =  1128  ft/s: 


f(Hz) 

Wavelenath 

20 

56.4  ft 

100 

11.3  ft 

500 

27.1  in. 

1000 

13.5  in. 

5000 

2.7  in. 
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7.  Direct  conversions; 


dB  =  1 0  log  (l/E-1 2)  [I  in  W/m**2] 
dB  =  20  log  (P/20E-6)  [P  in  Pascals,  RMS] 
dB  =  20  log  (P/29E-6)  [P  in  Pascals,  peak] 

P(Pascals  RMS)  =  20  E-6  *  antilog(dB/20)  =  20  E-6  *(10**(dB/20)) 

P(Pascals  peak)  =  28.9E-6*  antilog(dB/20) 

IIW/m**2]  =  ((Pascals,RMS)**2)/400 

l[W/m**2]  =  (Pascals, peak)**2/(2*D*V)  =  (Pascals,peak)**2/825.6  (air) 

1  Pascal  =  Newton/m**2  =  10  dyne/cm**2 
1  psi  =  6895  Pascals 

1  joule  =  E7  ergs  =  E7  dyne-cm  =Newton-meter 
1  atm  =  1.01325E5  Pascals  =  14.7  psi=1.013E6  microbars  (dyne/cm^) 

1  bar  =  E5  Pascals 
1  Watt  =  joule/sec 
1  Newton  =  E5  dynes 

8.  Technical  notes: 

Intensity  (RMS)  =[Pressure(RMS)]**2/(D*V) 

D=  air  density  (kg/m**3) 

V=  sound  velocity  in  air  (m/s) 

D*V  =  impedance  (at  20°C  and  1  atmosphere): 

D*V  (for  air)  =  1.2  [kg/m**3]  *  344  [m/sec]  =  412.8  [kg/(m**2-sec)] 

D*V  (for  water)  =  1000  *  1480  =  1,480,000  [kg/(m**2-sec)] 

For  the  same  acoustic  pressure,  the  intensity  in  air  is  (1 .48E6/41 3)  =  3,580  greater  than 
in  water. 

For  the  same  frequency  and  particle  displacement,  the  intensity  in  water  is  3,580  times 
greater  than  in  air. 

The  displacement  amplitude  of  air  =  A  =P(max)*V/(2*PrB*F). 

P(max)  =  maximum  (not  RMS)  pressure  amplitude.  Use  30  Pascals. 

B  =  air  adiabatic  bulk  modulus  at  normal  conditions  =  Gamma*air  pressure  =  1.4* 
1.013E5  Pascals  =  (142,000  Pascals). 

PI  =  3.14159  and  use  F=1 000  Hz. 

Gives:  A  =  1.2  E-5  meters  =  0.012  mm  (extremely  small). 

Wavelet  transforms  is  an  alternate  method  to  the  FFT  for  obtaining  good  frequency  resolution 
(long  recording  time)  up  to  higher  frequencies  (short  sampling  time)  for  pulses  and  portions  thereof. 
Here,  a  signal  region  is  expanded  on  an  ‘analyzing  wavelet’  function  that  is  relevant  to  the  signal 
under  consideration.  Examples  of  these  transforms  are  not  included  in  this  report. 
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AURAL  LIMITS  BY  MIL-STD-1474C 

A  distinction  is  made  between  steady  state  and  impulse  noise.  If  a  steady  state  noise  is  below 
85  dB  or  impulse  noise  is  below  140  dB  (determined  from  the  peak  pressure),  no  protection  is 
needed. 

Impulse  noise  is  characterized  by  two  parameters:  the  peak  pressure  level  in  dB  (dB  =  20*log 
(Peak  Pascal)/  29E-6)  and  the  effective  time  duration  (B-duration)  of  the  sound  in  milliseconds.  The 
dB  and  B  duration  lie  in  some  area  of  figure  12  and  by  table  4,  three  of  these  areas:  X,  Y,  and  Z 
require  specific  hearing  protection. 


Table  4 

Maximum  daily  exposure  to  impulse  noise 


Region  No  protection  Either  plugs  or  muffs  Both  plugs  and  muffs 


W  - Unlimited  Exposure - 

X  0  2,000  40,000 

Y  0  100  2,000 

Z  0  5  100 


A  single  exposure  consists  of  either  a  single  impulse  for  non-repetitive  system,  or  a  burst  as 
for  an  automatic  weapon. 

The  sequence  of  steps  to  determine  and  A  and  B  durations  will  be  illustrated  with  the  impulse 
time  trace  of  figure  13,  which  is  the  same  as  figure  2.  The  20  dB  intervals  are  marked. 

1 .  Find  the  peak  pressure  (voltage)  of  the  time  signal.  Peak  =  0.215  V. 

2.  Find  A,  time  for  pressure  to  go  from  zero,  rise  to  first  principal  positive  peak,  and  return 
to  zero.  A=  546  ms -536  ms  =  10  ms. 

3.  Set  two  bands,  L+  and  L-,  at  20  dB  below  the  peak  (10%  of  peak),  and  its  mirror  image 
about  zero  level.  L+  =0.021 5,  L-  =  -0.0215  V. 

4.  Find  B,  time  from  first  zero,  to  the  time  when  pressure  finally  stays  between  L+  and  L-. 
B  =  633  ms  -  536  ms  =  97  ms.  Q  =  extreme  time  point  633  ms. 

5.  Find  T 1 .  First  estimate  is  T 1  =5*A  =  5*1 0  =  50  ms. 

6.  However,  T1  must  be  equal  or  less  than  30  ms.  And  also,  if  T1  >  B,  then  start  reducing 
the  integer  from  5  to  4,3,2... 

7.  Sequence:  T1=4*A  =  40ms  (Still  greater  than  30  ms). 

8.  Sequence:  T1=3*A  =  30ms  (At  least  equal  to  30  ms.  Stop.) 

9.  Further  test:  Find  0.3*B.  0.3*B  =  0.3*97  =  29  ms. 

1 0.  Now  if  T1  is  equal  or  greater  than  0.3B,  then  B  =  duration  of  primary  portion. 

1 1 .  But  if  T1  is  less  than  0.3B,  find  the  first  point  past  0.3B  during  which  pressure  lies 

entirely  between  L+  and  L-.  Let  this  time  be  Q.  B  duration  is  then  Q-536ms.  In  this 
case,  T1  >  0.3B. 

12.  If  there  are  any  subsequent  fluctuations  that  vary  outside  the  limits  of  L+  and  L-,  and 
whose  time  duration  is  greater  than  10%  of  B,  then  add  this  time  duration  onto  B.  For 
this  trace,  there  are  none  and  the  B  duration  is  then  97  ms. 
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The  transformation  of  peak  volt  to  Pascals  was  0.21 5*1 2580.=2700  Pascal(peak).  Then 
dB=20*log(2700/28.9E-6)  =  159.4  dB.  By  figure  12,  this  dB  and  B  duration  fall  into  region  X. 


BEAT  FREQUENCY  FROM  MULTIPLE  ACOUSTIC  SOURCES 

The  usual  approach  of  considering  only  two  speakers  is  expanded  into  a  more  general 
arrangement  of  many  acoustic  sources  on  a  panel  perpendicular  to  the  ground,  and  having  micro¬ 
phones  located  anywhere  in  a  region  in  front  of  this  panel  array.  Some  calibration  must  first  be 
made  for  each  speaker  at  various  distances  to  obtain  at  least  a  geometric  regression  fit  for  sound 
received  at  each  mike.  The  input  data  needed,  and  results,  are  calculated  by  the  beat  frequency 
FORTRAN  program  in  table  5. 

The  acoustic  sources  are  considered  located  on  an  X,Z  plane,  perpendicular  to  the  ground, 
and  the  microphones  are  located  anywhere  at  XM,YM,  and  ZM,  where  YM  is  the  mike  ground 
distance  from  the  sources. 

The  program  assumes  there  is  complete  mixing  of  the  frequencies,  though  the  degree  of 
mixing  is  determined  experimentally.  The  frequency  content  at  the  location  of  a  microphone  in  the 
field  is  calculated  by  the  program.  This  content  is  independent  of  the  location,  though  there  is  a 
phase  difference  between  various  physical  locations. 

The  following  definitions  are  needed  to  describe  the  system  and  input  to  the  program: 

Xi=  X  coordinate  of  the  ith  speaker  in  X,M  plane; 

i=1 ,2,3,4... up  to  as  many  speakers  as  are  present. 

Zi=  2  coordinate  of  the  ith  speaker 

fi  =  frequency  of  the  ith  speaker  (frequency  can  be  any  explicit  function  of  time) 

XMj=  X  coordinate  (lateral)  of  the  jth  mike  in  the  vertical  XM,ZM  plane 
YMj=  y  coordinate  (axial)  of  the  jth  mike  from  the  plane  of  speakers. 

ZMj=  z  coordinate  (vertical)  of  the  Jth  mike  in  the  vertical  XM.ZM  plane 

C  =  speed  of  sound  =  1128  fps  at  20  C. 

Rij  =  actual  straight  line  distance  from  ith  speaker  to  the  jth  mike, 
t  =  time 

Pl=3. 141 592654 
i=  square  root  of  minus  one 

The  acoustic  signal  that  arrives  at  the  mike  has  a  peak  pressure  Si,  and  is  described  by: 

Signal  at  jth  mike  from  the  ith  speaker  =  Si  *exp[i*2*PI*fi*(t+Rij/C)] 

For  absolute  pressure,  the  mike  must  be  calibrated  with  the  speakers,  operating  at  some 
constant  frequency.  For  a  generic  expression  for  pressure  as  a  function  of  distance,  each  speaker 
is  calibrated  at  two  or  more  locations.  The  programs  uses  the  regression  coefficients  and  calculates 
the  signal  strength.  Si,  at  the  jth  mike  from  the  ith  speaker:  Si  =  A  *  (Rij  **  B)  where  A  and  B  are 
obtain  from  the  curve  fit.  If  only  two  calibration  points  are  obtained,  as  15  Pascals  at  a  distance  of 
5  ft,  and  2  Pascals  at  30  ft,  then  B  =  ln(2/15)  /  ln(30/5)  =  -2.01/(1.79)  =  -1.14,  and  A  =  15/(5**B)  = 
94. 
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The  frequency  content  is  added  to  the  signal  strength,  and  the  signal  is  now  expressed  in 
cosines  and  sines  from  the  complex  form: 

Signal  =  Si[  cos(2*PI*fi*(t+Rij/C))  +  i  *  sin(2*PI*fi*(t+Rij/C))] 

These  are  essentially  vectors,  and  can  be  added  as  vectors.  Figure  14  shows  the  addition  of 
three  such  sources  and  the  method  is  expanded  in  table  5. 

The  amplitude  of  the  resultant  vector  is  Ssum,  and  the  phase  angle  is  the  inverse  tangent  of 
the  ratio  of  the  imaginary  and  real  components.  The  straight-line  distance  between  any  speaker  j 
and  the  mike  i  is; 

Rij  =  SQRT(Q*Q  +  (DELX*DELX)) 

where  Q  =  SQRT(YMij*YMij  +  (DELZ*DELZ)) 

DELX  =  XMij  -  Xi 
DELZ  =  ZMij  -  Zi 


Table  5  is  the  FORTRAN  listing  of  the  program  illustrating  three  sources: 
f1=2000  Hz 

f2=  2100  Hz.  If  this  frequency  were  amplitude  modulated  at  30  Hz,  it  would  be  express 
ed  as  f2=  2100*  sin(PI*60*t),  but  such  an  expression  must  be  recompiled  in  the 
program. 
f3=2150  Hz 

For  these  three  acoustic  sources,  and  for  complete  signal  mixing,  the  following  figures  docu¬ 
ment  the  results: 

Figure  15  is  the  composite  time  signal  at  mike  1  for  two  source  signals;  2000  and 

2100  Hz. 

Figure  16  is  the  FFT  of  figure  15,  which  extracts  the  beat  of  100  Hz  and  harmonics. 
Figure  17  is  the  composite  time  signal  at  mike  1  for  three  source  signals;  2000,  2100, 

and  2150. 

Figure  18  is  the  FFT  of  figure  17,  which  extracts  the  beat  of  50,  100,  and  150  Hz  and 

harmonics. 

An  experimental  beat  frequency  from  2000  and  2100  Hz  sources  is  given  in  figure  19.  The 
total  recording  time  is  T=2  sec,  and  sampling  time  is  DT=20.833  ps.  The  frequency  range  is 
1/(2*DT)  =  24KHz.  Instead  of  choosing  the  96000  data  points  (T/DT),  for  a  frequency  resolution  of 
0.5  Hz  (DF=1/T),  the  VuPoint3  program  chose  2**17  data  points  (131072),  for  a  frequency  resolu¬ 
tion  DF=0.366  Hz. 


15 


The  FFT  of  the  positive  signal  is  shown  in  figure  20.  Whether  a  Hanning  or  no  filter  is  used  is 
irrelevant,  except  for  magnitude.  The  sum  and  difference  frequencies  are  extracted.  Figure  21  is  a 
finer  detail  of  the  carrier  frequencies  in  the  FFT.  Figure  22  is  the  integral  of  the  auto-correlation 
function,  graphing  the  cumulative  energy  of  the  various  frequencies.  The  zero  frequency  level  is 
approximately  32%  of  the  final  sum  in  figure  22.  When  this  zero  level  is  subtracted  out  of  the 
cumulative  sum,  the  percent  energy  contained  in  the  various  frequencies  becomes: 

Frequency  (Hz)  Energy  f%) 


100 

200 

2000 

2100 

4025 

4100 

4210 


11 

1 

37 

38 
1 

11 

1  Total=100% 


MODAL  FREQUENCIES  IN  AN  ENCLOSURE 

The  standard  expression  for  the  discrete  frequencies,  and  the  number  of  each,  has  been 
reduced  to  a  FORTRAN  program  in  table  6  for  a  rectangular  closed  structure.  For  a  closed 
reflecting  room,  the  fixed  wavelengths  that  can  exist  are  basically  given  as 

W=2/[SQRT((Nx/Lx)**2  +  (Ny/Ly)**2  +  (Nz/Lz)**2)] 

where  the  program  varies  the  integers  Nx,Ny,Nz  independently  from  0  to  38.  Lx,Ly,L.z  are  the  room 
dimensions.  The  longest  wavelengths  go  as  twice  the  room  dimensions,  and  for  one  open  end  they 
go  to  four  times  the  dimensions.  Figure  23  plots  the  specific  frequencies  that  can  exist  in  a  6.6  x 
9.3  X  20.6-ft  room,  and  the  number  of  each  up  to  a  frequency  of  220  Hz. 

To  illustrate  the  capability  of  finite  element  determination  of  acoustic  properties  in  an 
enclosure,  a  room  with  the  dimensions  described  in  figure  23  is  driven  by  an  acoustic  source  in  one 
corner  of  the  room  at  75  Hz  with  a  peak  to  peak  signal  of  two  Pascals.  Figure  24  plots  the  pressure 
contour  on  the  viewed  walls  at  one  instant  in  time.  The  resulting  pattern  is  not  simple. 

The  dB  pattern  of  this  pressure  is  plotted  in  figure  25.  To  illustrate  the  pressure  along  a 
diagonal  extending  from  the  source  to  the  other  end  of  the  room,  figure  26  shows  three  locations 
where  pressure  nulls  exist.  This  illustration  is  for  a  closed  room  with  reflecting  walls.  Variations  in 
closure  and  wall  absorption  are  readily  treated. 
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FFT.  Hanning.  Detrend.  KRUNlB.Cd 

FFTof Auto-Correlation  oF  Set  2  Uolts/hx 

(Uolts/hz]^2  N  VI  -0 


FFT.CmmEL  4.  PT=28.833iis  C.fl.G.  28  Oet  1997 
Cursor^  SetZ:  FFT.  Haimlng.  Detrend.  XRUN1B.C4;  Uolts 
Fnt«137  5.49983EH1  8  199.21907te 


Figure  3 

Fast  Fourier  transform  of  time  signal  pulse 


ENERGy.  FFT  AUTO.  CHAHNEL  4.  DT=28.833iis  C.H.G. 
Cursors  Set3:  FFTofftuto-Correlatlon  of  Set  Z:  Uolts 
PntSl37  3.BZ4B1E-7  B  199.21907hr 


Figure  4 

Energy  spectral  content  of  signal 
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INTEGRAL  OF  DfERGY.  FFT  MJTO.  CHfINNEL  4.  ST=ZB.833U 
Cursors  Set4:  Intg  of  FFTofAuto-Correlat:  Uolts<4iK 


Figure  5 

Integral  of  energy  content  of  signal 


INUERSE  OF  FFT.  CHANNEL  4.  DT=ZB.833us  Hanning. 
Cursors  Sets:  Inu_FFT  of  Set  Z:  Uolts 


Figure  6 

Recovered  time  signal  by  inverse  FFT 
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Figure  8 

Inverse  FFT  of  pulse  signal  without  time  extension 


KBUHIB.  CHAN. 4  COtIPLETE  TIHE  PULSE  8  TO  2  SET 

Cursors  Sell:  KnJMBie.bln  binary  data:  Uolts 
Pnt«48aei  -8.0393677  8  999.9984as 


Figure  9 

Acoustic  pulse  with  long  time  duration 


FFT.KRUNia.  CHMI.4  COMPLETE  TIME  PULSE  0  TO  2  SEC 
Cursors  Set8:  FFT  OF  ENTIRE  TIME  PULSE  0-2  SEC:  Uolts 
PntS547  6.94849E-4  0  199.95149hx 


Figure  10 

FFT  of  pulse  with  long  time  duration 
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C:SBINNSM1P1 

UOLTS 


zaee^ziea.zisa  hz.  »i=z.5e-s,  H=i8ee0.  6  hou  1997 
Cursors  Setl:  Col.Z:  C:SBINNSAnri:  VOLTS 
Pnt«5881  8.87955  8  lZ5.aiis 


SB.Bk 


lee.Bos  158.8k  zea.eus 

TIME 


Figure  17 

Composite  time  signal  for  2000,  2100,  and  2150  Hz 


B.aZ5 


FTI.  BEAT  Z888.Z108.Z1SB  HZ.  DT=Z.5E-5,  11=18888.  H 
Cursors  Set3:  FFT  OF  ZBes.  Z188  HZ.  HANHIMG.:  VOLTS 
PntSl  8.115766  8  8  te 


188. Bhz  Z8a.8hz  388.81iz 

PREQUENCy 


488.8112 


Figure  18 

FFT  of  composite  signals  of  2000,  2100,  and  2150  Hz 


ABEBZe.  FFT  T=2SEC. 07=20. 83US 
Uolts/hz 


ABE020.  T=2SEC.  BT==20.833US.  5  NOU  97.  €2^0,2100  HZ 
Cursor^  Set9:  0BE020.biii  binary  data;  Uolts 
Pnt«48627  0.0411987  0  1. 0130255s 


Intg  OF  AUTO  Correlation-ENEROy 

IUolt*«li*/|izr2  ABE020.  FFT  T=2SEC.DT=20.B3US 

N  >1  Uoli«/hz 


Figure  21 


Detail  of  2000  and  2100  Hz  signal  in  FFT 

MEB28.  CHftHWEL  1.  MAmiNG(«).  OMJLATIUE  EMERGY. 

OF"  ftUIO  Correlation-ENERGy;  Uolts«4iz 


Figure  22 

Cumulative  energy  versus  frequency  of  experimental  beat  signal 
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ANSYS  5.3 

FEB  13  1998 

13:30:43 

NODAL  SOLUTION 

STEP=1 

SUB  =1 

FREQ=75 

PRES 

SMN  =-.932052 
SMX  =1 


NFOR 


RFOR 


-.932052 

-.71738 

-.502707 

-.288035 

-.073362 

.14131 

.355983 

.570655 

.785328 

1 


f 


Figure  24 

Finite  element  determined  pressure  modes  in  room  driven  at  75  Hz 
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ANSYS  5.3 
FEB  13  1998 
13:56:59 

AVG  ELEMENT  SOLUTION 

STEP=1 

SUB  =1 

FREQ=75 

DB  (AVG) 

SMN  =44.074 
SMX  =90.39 


NFOR 


RFOR 


44.074 

49.22 

54.367 

59.513 

64.659 

69.805 

74.952 

80.098 

85.244 

90.39 


Figure  25 

Decibel  plot  of  pressure  modes  in  room  driven  at  75  Hz 
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Table  1 

Morse  complex  impedance  FORTRAN  program 


PROGRAM  MORSE 

C  ***DETERMINE  COMPLEX  IMPEDANCE  OF  MATERIALS**** 

C  ***PHILIP  MORSE,  'VIBRATION  AND  SOUND,'  1948,  PP. 239-246 
C  ***USE  COMPLEX  FUNCTION  OF  NATURAL  EXPONENT  TO 
C  ***OBTAIN  REAL  AND  IMAGINARY  PARTS  OF  IMPEDANCE.  10/24/97 
C  ******************************************************* 

C  ^ter  5  R^  data  on  one  line  in  file  MORSE  separated  by  commas. 
C  Have  the  file  MORSE  in  the  directory^  and  run  MORSE .EXE 
C  to  output  the  REAL  and  IMAGINARY  parts  of  the  impedance, 

C  and  the  ordinary  absorption  coefficient. 

c - 

C  1.  XLAM  -  wavelength  in  ft=sound  speed/frequency:  1128/50=22.56' 
C  2.  XL  -  distance  in  ft  between  speaker  and  specimen  face. 

C  NOTE:  Speaker  at  x=0;  specimen  face  at  x=L. 

C  3.  XA  -  X  coordinate  of  first  min  pressure  (from  specimen  face) 

C  4.  PMAX  -  maximum  pressure  (pascals  or  volts)  found  after  movino 
C  from  XA  toward  the  speaker.  ^ 

C  5.  PMINl  -  minimum  pressure  (pascals  or  volts)  at  location  XA. 

C  The  program  will  calculate  PI*ALPHA  and  PI*BETA  and  insert  them 
C  as  the  real  and  imaginary  components  of  S: 

C.  .S=PI*ALPHA  -  i*PI*BETA 

C..  FIND  ALPHA  AND  BETA  PROM  IMPEDANCE  TUBE: 

- - 

C.  .ALPHA  PROM:  TANH(PI*ALPHA)  =PRESS( MIN) /PRESS (MAX) 

C..BETA  FROM:  BETA=  1.-  2 . * ( XL-XA ) /XLAM 
C  XSI=DIMENSIONLESS  ACOUSTIC  IMPEDANCE  OF  MATERIAL 
C  XSI=  THETA  -  i*CHI  =  TANK ( PI*ALPHA-  i*PI*BETA) 

C  THETA  IS  THE  DIMENSIONLESS  MATERIAL  ACOUSTIC  RESISTANCE. 

C  CHI  IS  THE  DIMENSIONLESS  MATERIAL  ACOUSTIC  REACTANCE. 

C  THETA  =  l./MU,  WHERE  MU  IS  USED  IN  ANSYS  PROGRAM. 

C  ...  TAKE  AS  AN  EXAMPLE: 

C  MAX  PRESSURES  OCCUR  HALF  WAVELENGTH  APART. 

C  MIN  PRESSURES  OCCUR  HALF  WAVELENGTH  APART  ALSO. 

C  THEN,  DISTANCE  BETWEEN  MAX  AND  MIN  IS  ABOUT  QUARTER  WAVELENGTH. 

C .  NUMERICAL  EXAMPLE . 

C  XA  =  26  FEET  AS  LOCATION  OF  FIRST  MIN 
C  XLAM  =  22.56  FT  = (C/P=1128/50 ) . 

C  BETA  =  1-  2.* (XL-XA) /XLAM  =.64539 
C  PI*BETA  =  2.02753 

C  PI*ALPHA=  INVERSE  TANK  ( PMINl /PMAX ) =  ATANH( 0 . 1/0 . 2 ) 

C  PI*ALPHA=  0.5*LN( (1+PP)/(1-PP) )  ‘ 

C  WHERE  PP  =  PMINl/PMAX  =  0.1/0. 2  =0.5 
C  PI*ALPHA  =  .54931 
C  XLAM=22.56 
C  XL=30. 

C  XA=26. 

C  PMAX=.2 
C  PMINl =.l 
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Preceding  Page  Blank 


Table  1 
(cont) 


COMPLEX  S,XSI 
REWIND  1 

OPEN ( 1 , PILE* ' MORSE ' ) 

READ( 1,111)  XLAM,XL,XA,PMAX,PMIN1 
111  FORMAT (5F1 0.5) 

WRITE(*,222)  XLAM,XL,XA,PMAX,PMIN1 
222  F0RMAT(1X,'XLAM=',P9.4,'  XL=',F9.4,'  XA-',P9.4, 

1  '  PMAX=',P9.4,'  PMIN1*',P9.4) 

RATIO*PMAX/PMINl 

AC=1 .-( (RATIO-1 .)/ (RATIO+1 .)) **2 
PI=3. 14159265 


PIBETA*  PI*(1.-2.*(XL-XA) /XLAM ) 

C  PIALPHA  IS  INVERSE  HYPERBOLIC  TANGENT  OP  (PMINl/PMAX) 
PP»PMIN1/PMAX 


PIALPHA=0 . 5*LOG( ( 1 . +PP ) / ( 1 . -PP ) ) 

C  LOG  IS  TO  THE  BASE  E. 

S=CMPLX ( PIALPHA, -PIBETA ) 

C  FORTRAN  HAS  COMPLEX  EXPONENTIAL  FUNCTION  FOR  TANH. 
C  TANH*[E(S)  -  E(-S)]  /  [E(S)  +  E(-S)] 

C  E(S)  IS  E  TO  THE  POWER  S. 


100 

101 

102 

103 

104 


XSI=(CEXP(S)-CEXP(-S) )/(CEXP(S)+CEXP(-S) ) 

WRITE (*,100)  XSI 
WRITE (*,101) 

WRITE(*,102) 

WRITE(*,103) 

WRITE(*,104)  AC 

FORMAT (IX, 'THETA  ='  F10.4,'  CHI  -  P10.4) 

FORMAT ( IX, 'THETA=DIMENSIONLESS  REAL  PART  OP  IMPEDANCE.') 
FORMAT (IX,'  CHI=DIMENSIONLESS  IMAGINARY  PART  OP  IMPEDANCE.') 
FORMAT ( IX, 'THETA»ACOUSTIC  RESISTANCE,  CHI*ACOUSTIC  REACTANCE. 
FORMAT ( IX , 'ABSORPTION  COEFFICIENT  =  ' , P8 . 4 ) 

STOP 

END 


) 


Input  file  for  ‘MORSE’  or  ‘SWEN’: 

22.56^^26.^2^1 

Output  file  fi:om  ‘MORSE’  or  ‘SWEN’: 


XLAM*  22.5600  XL*  30.0000  XA-  26.0000  PMAX-  .2000  PMINl*  .1000 
THETA  =  1.2630  CHI  =  .7499 

THETA=DIMENSIONLESS  REAL  PART  OP  IMPEDANCE. 

CHI*DIMENSIONLESS  IMAGINARY  PART  OP  IMPEDANCE. 

THETA*ACOUSTIC  RESISTANCE,  CHI-ACOUSTIC  REACTANCE. 

ABSORPTION  COEFFICIENT  -  .8889 

Stop  -  Program  terminated. 
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Table  2 

Swenson  complex  impedance  FORTRAN  program 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  SWBN 

***DETERMINE  CXDMPLEX  IMPEDANCE  OP  MATERIALS 

SWENSON,  'PRINCIPLES  OP  MODERN  ACOUSTICS'  1965 
♦**USB  COMPLEX  FONCTIOH  OP  NATURAL  EXpSSl^O 

«=“  “a  imaginary  parts  of  the  lm|Sa^  “ 
and  the  ordinary  absorption  coefficient. 

1.  XIAM  “  in  ft=sound  speed/frequency:  1128/50=22.56' 

2.  XL  -  distance  in  ft  between  speaker  and  soecimen 

NOTE:  Speaker  at  x-0;  specimen  face  at  ^ 

I*  StAv  ^  at  first  min  pressure  from  the  specimen  face 

"frorL^SSarrSripe^ge?!""" 

*****»******i5**S?******f****i*»****°iJ>*®^^io®®jion  XA. 

Programs  cal^lates  PI*ALPHA  and  PI*BETA  and  inserts  them 
,  as  real  and  imaginary  components  of  S: 

C.  .S=PI*ALPHA  -  i*PI*BETA 

C. .PIND  ALPHA  AND  BETA  PROM  IMPEDANCE  TUBE. 

C  ****************************** 

TANH(PI*ALPHA)*  press  (MIN) /PRESS  (MAX) 

C..  BETA  PROM:  BETA=  1.-  ( (L-XA)/(XA-XB) )  ' 

C..XL=LENGTH  OP  IMPEDANCE  TUBE. 

^iRST  PRESS  (MIN)  PROM  SPECIMEN. 

COORDINATE  OP  FIRST  PRESS (MAX)  PROM  SPECIMEN  XB<XA 
C  XSI=DIMENSIONLESS  ACOUSTIC  IMPEDANCE  OF  MAtSiaL 
C  XSI*  THETA  -  i*CHI  ■  TANH(PI*ALPHA-i*PI*BETA) 

C  THETA  IS  THE  DIMENSIONLESS  MATERIAL  ACOUSTIC  RESISTANCE 
C  CHI  IS  THE  DIMENSIONLESS  MATERIAL  ACOUSTIC  REACTANCT 
C  THETA  =  l./MU,  WHERE  MU  IS  USED  IN  ANSYS  PRO^^ 

C  TAKE  AS  AN  EXAMPLE: 

C  L=30  FEET,  XA  =26. 

C  A=  (XA-L)/XLAM 

A^“?26-30)/22!sf’-''-'^?7?r™  “  50  HZ). 

TSaJ  ‘  *  5-36967 

COMPLEX  S,XSI 
REWIND  1 

OPEN ( 1 , PILE= ' SWEN ' ) 

READ(1,111)  XLAM,XL,XA,PMAX,PMIN1 
FORMAT (5P1 0.5) 

WRITE(*,222)  XLAM,XL,XA,PMAX,PMIN1 
F0RMAT(1X,'XLAM=',P9.4,'  XL=',P9.4,'  XA=',P9.4 
1  '  PMAX=',P9.4,'  PMIN1=',F9.4)  ' 

RATIO=PMAX/PMINl 

AC=1 .-( (RATIO-1 .)/ (RATIO+1 .)) **2 

XL=30. 

XA=26. 

XLAM=22.56 


111 

222 
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Table  2 
(cont) 


PMAX*0.2 

PMIN1=0.1 

A=(XA-XL)/XLAM 

PHI=3. 14159265* (1. -4. *A) 

C  CALCULATE  PHI,  AND  PLACE  THE  REAL  NUMBER  IN  SECOND  T  POSITION; 
S=CMPLX(0. ,-PHI) 

C  exponential  FUNCTION  FOR  TANK. 

C  TANH=[E(S)  -  E(-S)]  /  [E(S)  +  E{-S)l 
C  E{S)  IS  E  TO  THE  POWER  S.  '  ' 

RATIO*PMAX/PMINl 

Q» ( RATIO-1 .)/( RATIO+1 . ) 

XSI=(1.+Q*CEXP(S) )/(l.-Q*CEXP{S) ) 

WRITE(*,100)  XSI 
WRITE(*,101) 

WRITE(*,102) 

WRITE(*,103) 

WRITE (*,104)  .AC 

FORMAT (IX, 'THETA  ='  P10.4,  CHI  -  '  FIO  4 » 

IS?!^^m^^^'''™®^*“°^“^SIONLESS  REAL  PART  OP  IMPEDANCE  M 

IMAGINARY  PART  OP  IMPEDANCE  ') 
™ETA=AC0USTIC  RESISTANCE,  CHI  =ACOUSTIC  REACTANCE  '  1 
FORMAT (IX, 'ABSORPTION  COEFFICIENT  »',P8.4)  KiiACTANCE.  ) 

STOP 
END 


100 

101 

102 

103 

104 
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Table  3 

ASTM  complex  impedance  FORTRAN  program 


PROGRAM  ASTM 

C  ASTM  C384-77  IMPEDANCE  AND  ABSORPTION  OP  ACOUSTICAL 
C  MATERIALS  BY  IMPEDANCE  TUBE  METHOD.  (24  OCT  1997) 

Q  ****4r:ll;***********4r******4r***** ************************ 

C  Enter  7  REAL  data  on  one  line  in  file  ASTM  separated  by  commas. 

C  Have  file  ASTM  in  directory,  and  run  ASTM.EXE 
C  to  output  REAL  and  IMAGINARY  parts  of  impedance, 

C  and  the  ordinary  absorption  coefficient. 

C - 

C  1.  XLAM  -  wavelength  in  feet»sound  speed/f reguency :  1128/50=22.56' 

C  2.  XAC  -  X  coordinate  of  first  MIN  pressure  (from  end  plate  face) 

C  during  calibration  (metal  end  only  in  place;  NO  SPECIMEN) . 

C  3.  XCC  -  X  coordinate  of  second  MIN  pressure  (going  toward  speaker) 
C  during  calibration  (metal  end  only  in  place;  NO  SPECIMEN). 

C  4.  XA  -  X  coordinate  of  first  MIN  pressure  (from  specimen  face). 

C  5.  PMAX  -  MAX  pressure  (pascal  or  volt)  between  the  two  MINS. 

C  6.  PMINl  -  MIN  pressure  at  XA  (location  of  1st  min:  SPECIMEN) 

C  7.  PMIN2  -  MIN  pressure  at  XC  (location  of  2nd  min;  SPECIMEN) 

C  ***********  NUMERICAL  EXAMPLE  ******** 

C  STEP  1:  CALIBRATION.  METAL  PLATE  AT  CLOSURE  END  X=XL=30  FT. 

C  XA=FIRST  MIN  (NO  SPECIMEN)  -  XL-LAM/4  «  30-5.64  =  24.36 

C  XC=SECOND  MIN  (NO  SPECIMEN)*  XAA  -  XLAM/2  *  24.36-11.28  -  13.08 
C  XAC=24 . 36 

C  XCC=13.08 

C  XL* ( 3 . *XAC-XCC ) /2 . 

C  STEP  2:  SPECIMEN  PACE  IN  POSITION  AT  X=XL: 

C  REMEASURE  FIRST  AND  SECOND  MIN  PRESSURES:  (LOCATION  XC  NOT  NEEDED) 

C  XA*26. 

C  XC=14.72 

C  D1*XL-XA 

C  D2  IS  ABOUT  HALF  WAVELENGTH _  D2=XLAM/2  *  11.28 

C  D2=11.28 

C  PMINl, PMAX, PMIN2  ARE  PRESSURES  AT  XA,XB,XC: 

C  THESE  UNITS  CAN  BE  PASCALS  OR  VOLTS: 

C  PMAX*  .2,  PMINl  *  .1,  PMIN2-  .09 
REWIND  1 

OPEN ( 1 , FILE* ' ASTM ' ) 

READ( 1,110)  XLAM, XAC, XCC, XA, PMAX, PMINl, PMIN2 

110  FORMAT (7F8. 3) 

WRITE(*,111)  XLAM, XAC, XCC, XA 
WRITE (*,112)  PMAX , PMINl , PMIN2 

111  FORMAT (IX, 'XLAM, XAC, XCC, XA  *  ' ,5(1X,P8. 3) ) 

112  FORMAT (IX, 'PMAX, PMINl, PMIN2*  ' , 3 ( 1X,F8 . 3 ) ) 

XL* ( 3 . *XAC-XCC ) /2 . 

D1=XL-XA 
D2*XLAM/2 . 

C  USE  FORMULAS  11,13  IN  ASTM  FOR  ABSORPION  COEFFICIENT: 
RATI01*PMAX/PMIN1 
RATI02*PMAX/PMIN2 

ABl=l.-( (RATIOl-1. )/(RATI01+l. ) )**2 
AB2*l.-( (RATI02-1. )/(RATI02+l. ) )**2 
AB11=AB1- ( AB2-AB1 ) /2 . 

AB13=AB1-(D1/D2)*(AB2-AB1) 
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Table  3 
(cont) 


C  CONVERT  PRESSURE  UNITS  TO  DB: 

PMAX-20 . *LOG10 (PMAX/29E-6 ) 

PMINl =20. *LOGl 0 ( PMINl / 2  9 E- 6 ) 

PMIN2=20 . *LOG10 ( PMIN2/29E-6 ) 

XL1=PMAX-PMIN1 
XL2=PMAX-PMIN2 
XL0=XL1+ ( XL1-XL2 ) /2 . 

C  XKO  »  10**(XL0/20) 

XK0-EXP( (XL0/20)*ALOG(10. ) ) 

XM=(XK0+(1./XK0) )/2. 

XN-(XK0-(1./XK0) )/2. 

PHI=(D1/D2  -0.5) 

PHID=360.*PHI 

THETA«1 . / (XM-XN*COS ( PHI ) ) 

CHI*  ( XN*SIN  ( PHI ) )  /  ( XM-XN*COS  { PHI ) ) 

WRITE(*,100)  XL1,XL2,XL0,XK0 
WRITE(*,200)  D1,D2,XM,XN,PHI 
WRITE(*,300)  THETA, CHI 
WRITE (*,400)  AB11,AB13 

100  FORMAT(1X,'XL1,2,0,KO  , 4 ( 1X,F9 . 5 ) ) 

200  P0RMAT(1X,'D1,D2,M,N,PHI  =" ,5{1X,P9.5) ) 

300  FORMAT (IX, 'THETA,  CHI  =  ',  2F12.4) 

400  FORMAT(lX, 'ABSORPTION  COEFFICIENT  (EQ.ll ,13 ) »' , 2 ( 1X,P7 . 4 ) ) 
STOP 
END 


]iq}ut  file  for  ‘ASTM*: 


22.56, 24.36,13.(»^.^2^1r09 


Output  file  from  ‘  ASTM” : 


XLAM,XAC,XCC,XA  *  22.560  24.360 

PMAX, PMINl, PMIN2*  .200  .100 

XL1,2,0,K0  *  6.02060  6.93575 

D1,D2,M,N,PHI  *  4.00000  11.28000 

THETA,  CHI  =  1.8717  -.1858 

ABSORPTION  COEFFICIENT  (EQ.ll, 13)* 
Stop  -  Program  terminated. 


13.080  26.000 

.090 

5.56302  1.89737 

1.21221  .68516 

9053  .9005 


-.14539 
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Table  5 

Multiple  array  source  beat  frequency  FORTRAN  program 


PROGRAM  BEAT 

C  BEAT  FREQUENCY  FOR  MULTIPLE  POINT  ACOUSTIC  SOURCES  FOR  OPEN  REGION. 

C  FOR  RESPONSE  FOR  OP  TO  THREE  MICROPHONE  LOCATIONS. 

C  INPUT  FILE  IS  ARRANGED  AS  FOLLOWS ;( COMMA  SEPARATION  BETWEEN  DATA) 

C  1. NUMBER  OP  SPEAKER  SOURCES (MAX»20 ) ;  NUMBER  MIKES (MAX^S) 

C  2. THE  FREQUENCIES  OF  THE  SOURCES  (IN  ORDER) 

C  3.THE  X  COORDINATE  (LATERAL)  OP  THE  SOURCES  (ALL  IN  X,Z  PLANE). 

C  4. THE  2  COORDINATE  (VERTICAL)  OP  THE  SOURCES  (ALL  IN  X,Z  PLANE) 

C  5. THE  A  GEOMETRIC  REGRESSION  COEFFICIENT  FOR  THE  SOURCES. 

C  6. THE  B  GEOMETRIC  REGRESSION  COEFFICIENT  FOR  THE  SOURCES  (S-A*DIST**B) 
C  7. THE  XM  COORDINATE  (LATERAL)  OP  THE  MIKES. 

C  8. THE  YM  COORDINATE  (AXIS)  OP  THE  MIKES.  (SOURCES  ARE  AT  YM=0.) 

C  9. THE  2M  COORDINATE  (VERTICAL)  OF  THE  MIKES. 

C  10. NUMBER  OF  CYCLES  OF  FIRST  FREQUENCY,  AND  DATA  POINTS/CYCLE. 

C  11. SPEED  OF  SOUND. 

C++NOTE:  IF  SOUND  SPEED  IN  M/SEC,  ALL  DISTANCES  IN  METERS. (344  M/S) 
C++NOTE;  IF  SOUND  SPEED  IN  FT/SEC,  ALL  DISTANCES  IN  FEET.  (1128  FPS) 
C....IP  A  FREQUENCY  IS  A  FUNCTION  OP  TIME,  IT  MUST  BE  PROGRAMMED. 

C... EXAMPLE  OF  THE  SOURCE  PILE  WITH  THREE  SOURCES,  TWO  MIKES; 

C  3,2  : THREE  SOURCES,  TWO  MIKES 

C  2000,2100,2150  :FREQUENCIES  OF  THE  THREE  SOURCES. 

C  0,0,0  tLATERAL  X  COORDINATE  OP  SOURCES  1,2,3. 

C  0,4,8  ;VERTICAL  Z  COORDINATE  OP  SURCES  1,2,3. 

C  110. ,150. ,180.  :A  COEFFICIENT  FOR  SOURCES  1,2,3 
C  -1.13,-1.17,-1.22  ;B  EXPONENT  FOR  SOURCES  1,2,3 
C  0.,4.  ; LATERAL  XM  COORDINATE  FOR  MIKES  1,2 

C  100., 50.  :AXIS  YM  COORDINATE  FOR  MIKES  1,2 

C  2., 2.  ; VERTICAL  ZM  COORDINATE  FOR  MIKES  1,2 

C  500., 20.  ; NUMBER  CYCLES  OP  2000  HZ  WITH  20  DATA/CYCLE 

C  1128.  rSPEED  OP  SOUND  (HERE  1128  FT/SEC  FOR  20  C) 

C - PLOT  AMP1,AMP2,AMP3  FOR  RESULTANT  AMPLITUDE  AT  MIKES  1,2,3 

C - PLOT  PHAl , PHA2 , PHA3  FOR  RESULTANT  PHASE  (NOT  AS  SIGNIFICANT). 

DIMENSION  F(20) ,X(20),Z(20),A(20),B(20) ,XM(3) ,YM(3) ,ZM(3) 

DIMENSION  R(20),S(20),DELX(20),DELZ(20) 

REWIND  10 

OPEN ( 1 0 , PILE- ' BEAT ' ) 

READ(10,111)  NSOURCE,NMIKE 
READ(10,112)  (F(I) ,I=l,NSOURCE) 

READ(10,112)  (X(I),I*l,NSOURCE) 

READ(10,112)  (2(1) ,I*l,NSOURCE) 

READ(10,112)  (A(I) ,I-l,NSOORCE) 

READ(10,112)  (B(I) ,I=1,NSOORCE) 

READ(10,112) (XM(I) ,I-1,NMIKE) 

READ(10,112) (YM(I),I»1,NMIKE) 

READ(10,112) (ZM(I) ,I=1,NMIKE) 

READ(10,112)  CY,DCY 
READ(10,112)  C 

111  FORMAT (213) 

112  FORMAT (10F8. 3) 

113  FORMAT  (10F8. 2) 

WRITE (*,111)  NSOURCE , NMIKE 
WRITE(*,113)  (P(I) ,1=1, NSOURCE) 

WRITE(*,112)  (X(I) ,1=1, NSOURCE) 

WRITE(*,112)  (Z(I) ,1=1, NSOURCE) 
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Table  5 
(cont) 


WRITE (*,112) 
WRITE (*,112) 
WRITE (*,112) 
WRITE (*,112) 
WRITE (*,112) 
WRITE (*,112) 


(A(I),I=1,NS0URCE) 
(B(I) ,I*1,NS0URCE) 
(XM(I) ,I=1,NMIKE) 
(YM(I) ,I*1,NMIKE) 
(ZM(I),I=1,NMIKE) 
CY,DCy,C 


REWIND  1 
REWIND  2 
REWIND  3 
REWIND  4 
REWIND  7 
REWIND  8 

OPEN ( 1 , PILE- ' AMPl ' ) 

OPEN ( 2 , PILE- ' AMP2 ' ) 

OPEN ( 3 , PILE- ' AMP  3 ' ) 

OPEN( 4 ,PILE= ' PHAl ' ) 

OPEN ( 7 , PILE- ' PHA2 ' ) 

OPEN ( 8 , PILE- ' PHA3 ' ) 

PI-3.141592654 

C  NUMBER  OP  MICROPHONES  IN  FIELD 
DO  1  I-1,NMIKE 

C  CY-NUMBER  OP  CYCLES  OP  FIRST  FREQUENCY  P(l)  TO  BE  CONSIDERED  -500. 
C  DCY-POINTS  PER  CYCLE  FOR  DIGITIZING  -  20 
C. .DT-SAMPLING  TIME  -  1/(P(1)*DCY)  -  1/(2000*20)  -2.5E-5 
C. .  T-TOTAL  TIME  -  CY/F(1)  -  500/2000  »  .25  SEC 
C..  N-NUMBER  DATA  =  T/DT  -  .25/2.5E-5  -  10000 
N=CY*DCY 
DT=1./(P(1)*DCY) 

WRITE(*,400) 

400  FORMAT(15X, 'SIGNAL  STRENGTH  DISTANCE  R  BETWEEN',/ 

1  5X,'MIKE  SOURCE  S  AT  R  SOURCE  AND  MIKE') 

C  FIND  DISTANCE  OF  THE  SOURCES  AND  SOURCE  PEAK  TO  THIS  MIKE: 

DO  3  J=l,NSOORCE 
DELX(J)=XM(I)-X(J) 

DELZ(J)-ZM(I)-Z(J) 

RS-SQRT ( YM ( I ) * YM ( I ) +DELZ ( J ) *DELZ ( J ) ) 

R ( J ) -SQRT ( RS*RS+DELX ( J ) *DELX ( J ) ) 

C  CALCULATE  SIGNAL  SOURCE  AT  THE  MIKE  FROM  GEOMETRIC  REGRESSION  A,B: 
S(J)-A(I)*(R(J) )**B(I) 

WRITE(*,401)  I,J,S(J) ,R(J) 

401  FORMAT(2X,2I6,2(3X,F10.3) ) 

3  CONTINUE 

C  DO  TIME  CALCULATIONS 
DO  2  K=1,N 
T=DT*(K-1) 

C  SUM  UP  ALL  THE  CONTRIBUTIONS  OP  SIGNAL  SOURCES 
SUMREAL-0 . 

SUMIMAG-0. 

DO  22  KK=l,NSOURCE 

SREAL=S(KK)*COS(2.*PI*P(KK)*(T+R(KK)/C) ) 

SIMAG-S ( KK ) *SIN ( 2 . *PI*F ( KK ) * ( T+R ( KK ) /C ) ) 

SUMREAL-SUMREAL+SREAL 
SUMIMAG-SUMIMAG+SIMAG 
22  CONTINUE 
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Table  5 
(cont) 


C  FIND  AMPLITUDE  AND  PHASE  OF  COMPOSITE  SIGNAL  AT  MIKE  (I) 
IF(SUMIMAG.EQ.O. )  SUMIMAG-lE-8 
AMP*SQRT ( SUMREAL*SUMREAL+SUMIMAG*SUMIMAG) 

PHA= 1 8  0 . *2 . * ATAN ( SUMREAL/SUMIMAG ) / ( PI*  2 . ) 
IF(I.EQ.1)WRITE(1,500)  T,AMP 
IP(I.EQ.2)WRITE(2,500)  T,AMP 
IP(I.EQ. 3)WRITE( 3,500)  T,AMP 
IP(I.EQ.1)WRITE(4,500)  T,PHA 
IF(I.EQ.2)WRITE(7,500)  T,PHA 
IP(I.EQ.3)WRITE(8,500)  T,PHA 
500  FORMAT (2F10 .5) 

2  CONTINUE 

WRITE(*,402) 

402  FORMAT (IX, 50 ( IH* ) ) 

1  CONTINUE 

WRITE (*,403)  N,DT 

403  FORMAT ( IX, 'N,DT  = " ,16 , 1X,F12 . 6 ) 

STOP 

END 


Input  file  for  ‘BEAT’: 


Output  file  from  ‘BEAT*: 


3,2 

2000. ,2100. ,2150. 

0. ,0. ,0. 

0. ,4. ,8. 

110. . 150. .180. 
-1.13,-1.17,-1.22 

20. . 2. 

100. . 50. 

20. . 12. 

500. ,20. 

1128. 


.000  4.000  8.000 

110.000  150.000  180.000 
-1.130  -1.170  -1.220 

20.000  2.000 
100.000  50.000 

20.000  12.000 
500.000  20.0001128.000 

SIGNAL  STRENGTH 
SOURCE  S  AT  R 

1  .579 

2  .583 

MIKE  SOURCE  S  AT  R 
2  1  1.492 

2  2  1.519 


DISTANCE  R  BETWEE 
SOURCE  AND  MIKE 
103.923 
103.228 
102.684 

***************** 
DISTANCE  R  BETWEE 
SOURCE  AND  MIKE 
51.459 

2  3  -  50.675 

r^TTSoSo*******’*********************— >*•** 

«,DT  »  10000  .000025 

stop  -  Program  terminated. 
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Table  6 

Frequency  modes  in  closed  room  FORTRAN  program 


ON. 


] 


PROGRAM  MODES 

C  MODES  IN  SHOEBOX  ROOM.  D(X) ,D<y ) ,D(Z) 

C  ROOM  IN  SIZES  UP  TO  (X,Y,Z)  IN  INCREMENTS  IN  EACH. 

C  DO  FOR  SINGLE  SIZE:  X«9.32,  Y=6.562,  Z-20.6  FT 
C  SOUND  SPEED  C  TAKEN  AS  1128  PT/SEC. 

C  FIND  P(I,J,K)  FREQUENCIES,  AND  ASSOC  WAVELENGTH  W(I,J,K) 

C  W(X)*C/P(X)  W(I)«C/F(I) 

C  FOR  EACH  ROOM  CONDITION,  INTEGERS  NX,NY,NZ  VARY  FROM  1 
C  FOR  CLOSED  WALLS:  W(X)=2*D(X)/NX 

C  FOR  EACH  DIMENSION  D,  INTEGERS  GO  FROM  1,2,3,,,, 

C  CONDITIONS : 

C  W(I) .LE.2*D(I)/N(I)  OR;  W(X) .LE.2*D(X)/NX 

C  riF  ONE  END  CLOSED,  OTHER  OPEN:  W(I) .LE.4*D(I)/N(I) 

C  ...METHOD;  FOR  FIXED  SIZE  ROOM,  VARY  INTEGERS  FOR 
C  SIMPLE  FREQENCIES  IN  UNCLUTTERED  ROOM.. 

C  I,J,K  -  X,Y,Z  DIRECTIONS . DISTANCES  IN 

DIMENSION  DX(50) ,DY(50) ,DZ ( 50 ) ,NX( 50 ) ,NY( 50 ) ,NZ ( 50 ) , 

1  W(38,38,38),P(38,38,38) ,A(1000),B(1000) 

REWIND  1 
REWIND  2 
REWIND  3 

OPEN (1, FILE- 'FREQ' ) 

OPEN ( 2 , FILE- 'WAVE ' ) 

Qpgfj  /  3  FIIiE*  ^  FEEQN  ^  ) 

C  FILE  FREQ  GIVES  FREQUENCIES  ALLOWED,  AND  MODES  AT  EACH  FREQUENCY 
C  FILE  FREQN  GIVES  FREQUENCIES  ALLOWED,  AND  ACCUMULATED  NUMB^.  . 

C  FIND  NUtfflER  MODE  IN  2  HZ  INTERVALS,  UP  TO  IA*2  HERTZ  (lA  GROUPS) 
C  lA  MAXIMUM  »  1000 
IA-110 

DO  77  1=1, lA 
A(I)»0. 

B(I)-0. 

77  CONTINUE 

FMIN-100000. 

C-1128. 

C2-C/2. 

INDEX- 10 
DO  1  K=l,l 
C  DZ(K)=K 

DZ(K)-20.6 

KS=K 

DO  2  J»l,l 
C  DY(J)=J 

DY(J)=9.32 

JS=J 

DO  3  1=1,1 
C  DX(I)=I 

DX{I)=6.562 

IS=I 

C  VARY  THE  INTEGERS . . 

NUM-38 

DO  4  KK=1,NUM 
NZ(KK)=KK-1 
DO  5  JJ=1,NUM 
NY(JJ)=JJ-1 
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Table  6 
(cent) 


DO  6  II«1,NUM 
NX(II)»II-1 

FF=C2*(SQRT( (NX(II)/DX(I) )**2  +  (NY( JJ)/Dy ( J) )**2 
1  +  (NZ(KK)/DZ(K) )**2) ) 

IP(FP.GT.O.)  WW=2./(SQRT( (NX(II)/DX(I) )**2  +  (NY ( JJ)/DY( J) ) **2 

1  +  (NZ(KK)/DZ(K) )**2) ) 

C  STORE  AN  ACCEPTABLE  STANDING  WAVELENGTH 
W(II,JJ,KK)«WW 
F(II,JJ,KK)=FF 
200  PORMAT(P10.2) 

C  FIND  CONDITIONS  WHEN  LOWEST  FREQUENCY  OCCURS.. 

IF(FF.LT.FMIN)  GO  TO  50 
GO  TO  51 

50  XL«DX(I) 

YL*DY(J) 

ZL-DZ(K) 

IXN=II-1 

JYN=JJ-1 

KZN=KK-1 

FMIN»FF 

51  CONTINUE 

100  FORMAT (IX, 'ROOM,  K#,  WAVEL,  PREQ=  ' , 3I3,2X, 313 , 2 ( 1X,F8 . 2 ) ) 

WRITE(2,600)  WW,PP 
DELP*2 . 

DO  300  IN=1,IA 

DPB=DELF*IN 

DPA=DELF*(IN-1) 

C  A  GIVES  THE  FREQUENCY (MEAN  OF  1  HZ),  B  GIVES  NUMBER  MODES. 
IP(FF.GT.DFA.AND.PF.LE.DPB)  A(IN)»(DFB+DPA)/2. 
IF(FF.GT.DPA.AND.FF.LE.DPB)  B(IN)=B(IN)+1. 

300  CONTINUE 

10  CONTINUE 

6  CONTINUE 

5  CONTINUE 

4  CONTINUE 

3  CONTINUE 

2  CONTINUE 

1  CONTINUE 

BMODES«0 

9  DO  90  IN*1,IA 

C  DO  NOT  WRITE  IP  NO  FREQUENCIES  WERE  POUND; 

IF(B(IN) .EQ.O. )  GO  TO  90 
WRITE(1,600)  A(IN),B(IN) 

600  FORMAT (2F1 0.2) 

BMODES=BMODES+B ( IN ) 

WRITE(3,600)  A(IN),BMODES 
WRITE(*,400)  B(IN) ,BMODES,A(IN) 

400  FORMAT(lX, 'FREQ  INTERVAL,  NUMBER  IN,  TOT  MOD=  ' , 3(1X,F10.2) ) 

90  CONTINUE 

WRITE(*,222)  XL,YL,ZL,IXN,JYN,KZN 
222  FORMAT (IX, 'ROOM  SIZE;INDEX  FOR  LOWEST  PREQ=  ',  3F6 . 2, 2X, 314 ) 

STOP 
END 
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